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A CHARACTERIZATION OF ENERGY-PRESERVING METHODS 
AND THE GONSTRUCTION OF PARALLEL INTEGRATORS FOR 
HAMILTONIAN SYSTEMS 

YUTO MIYATAKE* AND JOHN C. BUTCHERt 

Abstract. High order energy-preserving methods for Hamiltonian systems are presented. For 
this aim, an energy-preserving condition of continuous stage Runge-Kutta methods is proved. Order 
conditions are simplified and parallelizable conditions are also given. The computational cost of our 
high order methods is comparable to that of the average vector field method of order two. 
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1. Introduction. This paper is concerned with the numerical integration of 
Hamiltonian systems of the form 

^y = fiy), y(to) = yoeK^, 

where f{y) = S'^H{'y) with a non-singular skew-symmetric constant matrix S and 
sufficiently differentiable Hamiltonian H : —>• R. Although the dimension N of the 

system is even, in the case of Hamiltonian systems, we do not assume this because the 
results are also applicable to Poisson systems with constant Poisson structure, where 
S is not necessarily non-singular and N need not be even. 

Several geometric properties are known in classical mechanics. For example, the 
exact flow is energy-preserving in the sense that 

AiJ(y) = = VH(y)^5ViI(y) = 0, 

as well as being symplectic. Numerical integrators exactly inheriting such properties 
are called geometric numerical integrators [11], and the most remarkable advantage 
of adopting geometric numerical integrators is that they often give qualitatively cor¬ 
rect numerical solutions over an extremely long period of time. Although integra¬ 
tors inheriting such geometric properties as much as possible would be preferable, it 
is known to be impossible to construct methods inheriting both symplecticity and 
energy-preservation [7, 23]. This is because if such methods exist the numerical flow 
coincides with the exact flow. In the last decades, methods inheriting one of these 
properties, i.e. symplectic methods and energy-preserving methods, have been devel¬ 
oped. 

Although both symplectic methods and energy-preserving methods have their own 
advantages, much more attention has been devoted to symplectic methods. One of the 
most remarkable advantages of symplectic methods is that they also nearly preserve 
the energy without any drift. It should also be noted that there is a symplectic 
condition for Runge-Kutta (RK) methods [15, 19, 20]: 
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where s is the stage number, and and bi are the coefficients of RK methods. This 
condition is not only sufficient but also necessary under a certain assumption. In fact, 
symplectic implicit RK methods exist whose typical examples are Gauss methods. 
For more details, we refer the reader to [11]. 

While symplectic methods are useful for most Hamiltonian systems, it is some¬ 
times mandatory to adopt energy-preserving methods. For example, errors in the 
Hamiltonian obtained by symplectic methods can sometimes be fatal for chaotic dy¬ 
namics. Thus, it is strongly hoped that the study on energy-preserving methods 
reaches to the same maturity as symplectic methods. However, the study on energy¬ 
preserving methods is more challenging and, in fact, newer than that on symplectic 
methods. One reason is that no RK method is energy-preserving for general Hamil¬ 
tonian systems [5]. 

The projection method is a relatively simple method of constructing energy¬ 
preserving integrators. The idea is very simple, but the implementation and the 
construction of the projection is specific to individual problems. An even more se¬ 
rious disadvantage of projection-based integrators is that the long term behaviour 
usually deteriorates. Compared with the projection method, the discrete gradient 
method, which was first proposed by Gonzalez [9] (see also McLachlan et al. [16]), is 
rather a more sophisticated and systematic approach. 

After some pioneering studies, substantial progress was made in the last decade. 
From a theoretical viewpoint, Faou-Hairer-Pham [8] and Chartier-Faou-Murua [7] re¬ 
vealed that there exist energy-preserving B-series integrators other than the exact flow 
by showing an energy-preserving condition. Furthermore, concrete energy-preserving 
methods have also been developed, Quispel-McLaren proposed the average vector 
field (AVF) method [18]. This method is energy-preserving, a subclass of the discrete 
gradient method, a B-series method and of order two. A more recent interest is to find 
a high order energy-preserving method as an extension of the AVF method. Hairer 
proposed the AVF collocation method [10], and Brugnano et al. proposed the Hamil¬ 
tonian boundary value method [1]. These methods are essentially the same methods, 
although the latter method has been developed mainly for polynomial Hamiltonian 
systems. These methods are based on so called continuous stage RK (CSRK) meth¬ 
ods. It should be noted that roughly speaking the computational cost of the AVF 
collocation method is almost the same as the same order Gauss method. 

The main aim of this paper is to construct a more efficient energy-preserving 
method than exists at present. To achieve this aim we focus our attention on a matrix 
M which characterizes GSRK methods. In terms of M, we consider (i) an energy 
preserving condition, (ii) order conditions and (iii) criteria for parallel implementation 
using real arithmetic. We then put together these three perspectives to construct high 
order energy-preserving methods which can be computed in parallel. We then derive 
concrete integrators of up to order six. 

The paper is organized as follows. In Section 2, some existing energy-preserving 
methods are briefly reviewed with the introduction of CSRK methods. In Section 3, 
an energy-preserving condition of CSRK methods is discussed. In addition to the 
sufficient condition already given in [17], we show that the condition is also neces¬ 
sary under a small assumption. The condition is expressed in terms of the matrix 
M. In Section 4, order conditions of energy-preserving CSRK methods are character¬ 
ized in terms of the matrix, where simplifying assumptions play an important role. In 
Section 5, parallel energy-preserving methods are constructed. We first consider a par- 
allelizable condition of CSRK methods, and then construct parallel energy-preserving 
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methods. Furthermore, we investigate fourth-order methods in detail. The derived 
fourth order methods can be implemented faster than the same order AVF collocation 
method, but an actual error of our methods per each time step becomes bigger than 
that of the existing method. We therefore have to discuss to what extent our method 
is efficient. Finally, concluding remarks are given in Section 6. 

2. Energy-preserving continuons stage Runge Kutta methods. In this 
section, we review some existing energy-preserving methods and continuous stage 
Runge-Kutta methods. 

The discrete gradient method proposed by Gonzalez [9] (see also McLachlan et 
al. [16]) is a first systematic method for deriving energy-preserving integrators. The 
key idea of the method is to introduce a sophisticated approximation of the gradient. 
We first define a discrete gradient, a map VR : x —)• K^, satisfying 

(2.1) H{x) - H{y) = VH{x, y)^{x - y), 

(2.2) yH{x,x) = H{x) 

for any x,y € , and then define an integrator by 

(2.3) yi^ = SVH{y,,yo), 

where h denotes the stepsize. Here, the condition (2.1) corresponding to the chain 
rule ^H{x) = ViJ(x)^x is called the discrete chain rule, and the condition (2.2) 
requires the consistency. Due to the discrete chain rule and skew-symmetry of the 
matrix 5, the energy-preservation for the discrete gradient method (2.3) can be easily 
verified: 

^{H{yi)-H{yo)) = VH{y,,yoV^^^=VH{y,,yoVSVH{y,,yo)=0. 

The discrete gradient is not generally unique and several approaches have been pro¬ 
posed. The average vector field (AVF) gradient proposed by Quispel-McLaren [18] 

VH{yi,yo)= f VH{ryo{1 - r)yi) dr 
Jo 

is one of several possible approaches. In general, the AVF method reads 

(2.4) yj^=yo + h [ /(rj/o + (1 - T)yi) dr. 

Jo 

The AVF method is a symmetric, second order and furthermore B-series method as 
will be shown later. 

After the proposition of the AVF method, the AVF method was extended to higher 
orders. Hairer extended the AVF method by slightly modifying the idea of the stan¬ 
dard collocation methods [10] (see also the Hamiltonian boundary value method [1]). 
In this paper, we call the extended method the AVF collocation method. As is the 
case with the standard collocation method which can be interpreted as a RK method, 
the AVF collocation method can be interpreted as a continuous stage Runge-Kutta 
(CSRK) method. 

Definition 2.1 (Continuous stage Runge-Kutta methods). Let he a poly¬ 
nomial in T and (). Assume that = 0. We denote by s the polynomial degree of 
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At in T. Let Bq be defined by = Ai We search for an s-degree polynomial Y^ 
(t G [0,1]^ and yi such that they satisfy 


(2.5) 


Yr =yo + h f Arxf{Yc) dC, 
^0 

yi=yo-\-h [ Brf{Yr) dr. 
^0 


A one-step method yo yi is called an s-degree continuous stage Runge-Kutta 
(CSRK) method. A CSRK method is said to be consistent if fg Bi^ d( = 1. 

Remark 1. CSRK methods were originally due to Butcher [2] (see also [3]), 
where A^x was not assumed to be a polynomial but to be bounded. In this paper, 
we focus on a polynomial Arx because this restriction simplifies discussion below and 
more general functions do not seem to be of greater benefit to our aim than polynomial 
functions. 

In the above definition, = 0 is assumed so that Yq = yo. Y^ can be seen as 
an approximation of y(to + Crh), where Cr = Jq dC- The final stage yi is also 
expressed by yi = Yi due to the relation Since in this paper we shall focus 

on which is polynomial of degree s in r and s — 1 in we often express A,-.^ as 


( 2 . 6 ) 


r^rX — 


M 


1 

c 

>s-l 


with a constant matrix M G 

In the case of the AVF collocation method, the polynomial A^,^ is given by 


where 


S 1 

^rx = ^-r i Li{a)daLi{C), 

i=l -^0 


= n 


r — Ci 


. , . / . Cf Cj 


bi= Li{T) dr 

Jo 


and Cl,..., Cg are distinct, real numbers. It was proved in [10] that the AVF collocation 
method is energy-preserving for Hamiltonian systems independently of the Ci values. 
If the Ci values are the zeros of the s-th shifted Legendre polynomial, then the method 
has order p = 2s. Here lists concrete expressions of Arx for s = 1, 2, 3: 


s = 1 : A^x = T, 

s = 2 : At,c = t((4 - 3t) - 6(1 - r)C), 

s = 3 : Arx = t{{9-18t + lOr^) - 12(3 - 8 t -h Sr^)^ -b 30(1 - 3 t -b 2t^X^) . 

Note that in the case of s = 1 and Arx = t, the corresponding CSRK method 
coincides with the AVF method (2.4). 

The following theorem ensures the unique solvability of CSRK methods. 
Theorem 2.2 (Unique solvability [21]). Suppose that f satisfies a Lipschitz 
condition 


\\f{v)-f{v)\\ <L\\r]-r]\\ 
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with a constant L. If 


h < 


L(^max^g[o,i]/p 


there exists a unique solution of (2.5). 

This theorem can be proved in a similar way to the proof of the unique solvability 
of RK methods (see, e.g. [3, 12]) as mentioned in [21]. For the convenience of the 
reader, this proof is included below. 

Proof. First, we note that the function space of s-degree polynomials with the 
property Yq = yo can be identified with We define a metric on by 


p{Yr,Yr) 


max IIW 

rG[0,l] 




and consider the map 


Yr (j){Yr) = yQ + h [ A^xfiYi;) dC. 

Jo 


Then we estimate p{cj){Yr), cffYr)) as follows: 


p{(j){Yr),(j){Yr)) = h max 
re[o.i] 


< h max 

re[0.1] 


[' ArxifiYc) - f{Yc)) dC 

Jo 

["lArxWlfiYc) - f{Yc)\\ dC 
Jo 


</iL max [ |At- r| IIW — Fr II dC 
< hLl max / |ylT-A|dC| max ||Fc“^cII 

V^e[o.i]7o J Celo.i]" 

= /iL(max [ |ylr.,^|dc)/9(Fr,Fr)- 

V'>-e[o.i]yo / 


If the stepsize h satisfies /iL^max,-g[o_i]|^T,c|dC) < li the contraction mapping 
theorem ensures the existence and uniqueness of the fixed point of the map (j). □ 

As pointed out in [10], the numerical solution of a CSRK method can be inter¬ 
preted as a B-series. Below we introduce some notation to be needed in the subsequent 
sections. See [3, 11] for more details on B-series. 

Let T be the set of rooted trees 




}• 


We denote by |t| the number of vertices of a tree t. A new tree obtained by connecting 
the roots ofti,... ,tm G T is denoted by t = [ti,..., tm]- If some oiti,... ,tm equal to 
each other, we write t = [ff-,... ,^5]^], e.g. [.,1,1] = [•,! ]. The symmetry coefficient 
cr : T ^ R is defined recursively by 

m 

cr {. t ) = rfc!(T(4)’'^ 

k=l 


cr(,) = I 
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where t = ... ,^5^]. The elementary differential is a mapping F{t) : —>■ 

defined recursively by 

F{.){y) = f{y), Fit)iy) = f^^HF{FKy),F{t^){y)). 

For a map a : T U {0} —5- R, a formal series of the form 

B{a, y) = a{%)y + ^ ^a{t)F{t){y) 
teT ^ 


is called a B-series. The exact time-/i flow of a differential equation y = f{y) can be 
expressed as y(to + h) = B{e,y{to)), where the coefficients e are given by 


e(0) = e(,) = 1, e{t) = |^e(ti) • • • e{tm)- 

As is the case with Runge-Kutta methods, every CSRK method is a B-series integrator 
yi = B{(j),yo), where the coefficients (j), called the elementary weights, are given by 
0(0) = 0(,) = 1, 4’t{») = Cr and 

4>r{t) = f ArX<l>c{tl) ■ ■ ■ 0(0= f 5t0t (^ 1 ) ’ ' ' 0t (^m) dr. 

Jo Jo 

Here are the first four elementary weights: 

0(.) = [ Brdr, 

Jo 

0(Y) = [ BrArxAr^K drdCdK, 

Jo 


0 ( 1 ) = [ BrArxdrdC, 

Jo 

0(i) = [ Hr drdCdK. 

Jo 


In addition to the set of trees 7“, we define a set of forests F and a set of free 
trees FT- The forest is an unordered finite collection of trees from F : 

F = '{c,**,!,•••,•!, , J,... 

An element of < € is written as t = tp’’ {ti,... ,tp & J~), and the sum of the 

vertices are denoted by |t|, e.g. |,I^| = I.HI = 5. For t gF, B-{t) denote the forest 

consisting of the subtrees of t, e.g. H_(*^) = ,^I. For a forest t = ■ ■ -tp”, we 

define a map a : —>■ R by a(t) = a{tiY^ ■ ■ ■ a{tpY^. 

Next, we introduce the set of free trees FF. Note that there is a group of trees 

whose elements have topologically the same shapes, e.g. the shapes of 


T 


and I are the same. We define the set FF by using the Butcher product, for two 
trees u = [mi, ..., Uq\ G F and v gF, 


uo V := [ui,... ,Uq,v] G F. 


Note that Mou 7 f:uouin general. A tree u o v can be obtained by shifting the 
root of u o u. Multiple shift of root induces an equivalence. All trees in the same 
equivalent class has the same number of vertices, and have topologically the same 
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shape. We call each equivalent tree a free tree and denote the set of free trees by 
TT- The canonical projection is denoted by tt : T TT- This map is not bijective, 




e.g. TT 1, I}. For two equivalent trees u and v, we define a 

distance k{u,v) by the minimum number of root shifts necessary to obtain u from v. 


(f.T, 


For example, k{u, u) = 0, k{u, v) = k(v, u) and k( J , I) = 3. Any free tree containing 
a member with a factorisation of the form m o it is called superfluous tree. 


3. Necessary and Sufficient condition for energy-preservation. In order 
to construct efficient energy-preserving CSRK methods, it is mandatory to charac¬ 
terize an energy-preserving condition for CSRK methods. The sufficient condition 
was already given in [17]. Although one can proceed to the subsequent sections by 
using the condition, we also show that the sufficient condition is also necessary un¬ 
der a small assumption. The necessity indicates that it is almost impossible to find 
energy-preserving CSRK methods outside the framework presented in this paper. 

3.1. Sufficient condition. Theorem 3.1 ([17]). A CSRK method is energy¬ 
preserving, if the matrix M in (2.6) is symmetric, i.e. if is symmetric. 

For readers’ convenience, we here show the proof of this theorem. 

Proof. If the matrix M is symmetric, it can be diagonalized by an orthogonal 
matrix P: M — P^AP, where A is a diagonal matrix whose elements are given by 
Aii = Xi. The partial derivative -^Ar^ is expressed as 

^ A*(P<p(r))^(Pv3(C))i, 


where (^{t) = [l,r,... ,t® ^]^ and {■)i denotes the i-th component. By using these 
notation, we have 

H{y,)-H{yo)= Ai/(y^)dT = ^ Y^VH{Yr)dT 

= ^^-^C'5Vi7(n)dc) VH{Yr)dT 

= (^J\ppiO)^^HiYc) dc) (^j\pipiT)),VHiYr) dr) 

= 0 . 


The last equality follows from the skew-symmetry of S. □ 

Remark 2. Tang and Sun have also provided essentially the same sufficient 
condition in a different manner [22]. 

3.2. Necessary condition. The idea of proving the necessity is to use a neces¬ 
sary and sufficient condition of energy-preserving B-series methods shown by Chartier 
et al. [7]. We will also use some lemmas which are motivated by Celledoni et al. [6]. 

The necessary and sufficient condition of energy-preserving B-series methods is 
summarized in the following theorem. 

Theorem 3.2 ([7]). A B-series B{a,-) is energy-preserving if and only if 

-— a(B-(u)) = 0, for all t € J-T, 


(3.1) 



where t is a member ofTr~^{t). 

As corollaries, we get the following properties. 

Lemma 3.3. If a consistent B-series method is energy-preserving, then it follows 

that 

(3.2) a([/-i])(^=^ i for alike n+{~ {1,2,i,...}). 


Proof. The idea of the proof is similar to that of Theorem 2.2 in [6]. This theorem 
is proved by considering a free tree such that k vertices are connected to a common 
vertex, e.g. for fc = 4. Any free tree of this type has two equivalent trees: ti = 
and t 2 = [[•^~^]]- The symmetry coefficients are easily calculated to be cr(ti) = k\ 
and a(t 2 ) = {k — 1)!. Note that a{B-(ti)) = a(,)^ = 1 and a{B-{t 2 )) = 
Substituting a and a into (3.1), we obtain the condition (3.2). □ 

Lemma 3.4. If a consistent CSRK method is energy-preserving, then it follows 

that 


(3.3)p / f BrCP-^ArxC^dTdC-q f [ BrC^-^ArxCPdrdC 
Jo Jo Jo Jo 

for all p,q e N'*'. 


1 


9+1 p+ 1’ 


Proof The idea of the proof is similar to that of Lemma 2.3 in [6]. To prove 
this condition, we introduce double bush trees. As illustrated in the picture below, a 
double bush tree tp^q is a free tree having p leaves on the left side and r leaves on the 
right side (this figure is an example with p = 2 and 9 = 3). 



It is clear that tp^p is superfluous, and tp^q = tq^p. 

There are four distinct but equivalent rooted trees for tp^q-. ti = 
t 2 = [•^[•‘^]], ^3 = [[•^]»‘^] and t 4 = [[[•^]»‘^“^]] (it is clear that K{ti,tj) = \i-j\). Here 
are examples for t 2 , 3 . 


'r*', Jf. 


Thus, it follows that 

a{B_{ti)) 
a{B-{t2)) 
a{B-{t3)) 
a(H_ (t4)) 


= a([.^>-i, [.«]])= r r BrCr^ArxCUrdC, 

Jo Jo 

= «(.)M[.^]) = ^, 

= a(.)M[."]) = ^, 

p + 1 

= a([.«-\ [.P]]) = r r BrCr^ArxCPdTdC. 

Jo Jo 
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The second and third are due to Theorem 3.3. The symmetry coefficients are calcu¬ 
lated to be 


cr{t 2 ) = aits) = p\q\, a{ti) = p\iq - 1)\. 


Substituting k, a and a into (3.1), we obtain the condition (3.3). □ 
We consider an s x oo matrix $ defined by 


= 


_i-l 


a dr. 


We now show the following theorem. 

Theorem 3.5. If a consistent CSRK method for which $ is of full rank s is 
energy-preserving, then M = . 

Proof. We see by contradiction that the consistency and the full rank assumption 
indicate Br = C).. Under the assumption Br = C)., we can rewrite the condition (3.3) 
as follows. Using the integration-by-parts formula, we have 



Similarly, 


Substituting these relations into the condition (3.3), we have 

r r C^A'^^^CUrdC = f r drdC 

Jo Jo Jo Jo 


= cp[i,' 


1]M 




Note that 
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and 



which is equivaient to 


= 0 . 


Since we assume that the matrix $ is of fuii rank s, it foiiows that M = . □ 

Remark 3. Let us discuss the assumption about the rank o/$. The technique 
of the proof of Theorem 3.5 is motivated by the proof of the necessary condition of 
symplectic RK methods, especially Theorem 7.10 in [11, Chapter VI], where a similar 
sxoo matrix characterizes the reducibility of RK methods. In contrast to RK methods, 
however, we believe that for any consistent CSRK methods the rank of the matrix $ 
is always of full rank s. This conjecture can be easily verified for monotonic Cr or the 
case the function space spanned by {C'*}jgpj+ is dense in {v \ v € C[0,1], t>(0) = 0}, 
but the conjecture for more general cases has not been proved. We leave the proof of 
this conjecture to our future work. 

4. Families of energy-preserving continuous stage Runge Kutta meth¬ 
ods and order conditions. In this section, we characterize order conditions of 
energy-preserving CSRK methods in terms of the matrix M (2.6), based on the sim- 
piifying assumptions 


B{p): 

[ BrC’[-^AT = 

Jo 

1 

fc’ 

fc = 1, . . 


Civ): 

[ ki.,cC'^^-idC = 
^0 

k ’ 

fc = 1, . . 

■ ,v, 

m)-- 

[ BrC^-^Ar.cdT = 

Jo 

^(1 - cf). 

fc = 1, . . 
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and construct several families of energy-preserving methods. 

The order of a CSRK method satisfying the simplifying assumption B{p)^ C{r]) 
and D{^) is at least min(/ 9 , 2?7 -|- 2 ,77 ^ 1). The following two theorems can further 

simplify the order conditions. 

Theorem 4.1. For any consistent CSRK methods with the property M = , 

the B{p) condition is satisfied for all p = 1,2^.... 

Proof. For any CSRK methods with the property M = it follows thatSr = 
Cfi The consistency indicates that Ci = B^. dr = 1. Hence, it follows that for any 

k 

□ 

Theorem 4.2. For any consistent CSRK methods with the property M = , 

the < 7 ( 77 - 1 - 1 ) condition is equivalent to the D{r]) condition. 

Proof. We prove 



is equivalent to 


(4.2) BrC^-^ArxdT=^{l-C^). 

Due to M = , it follows that <7). = Br. Recall that we always assume that 

(7o = /g dC = 0, and the consistency indicates (7i = 1. Let 


At.C, — 


M 


' C ■ 
C72 

CVs 


so that and Ar^c^ = Aq^t- 

We start with (4.1). By using the integration-by-parts formula, we see that (4.1) 
is equivalent to 


Ar.cCt 


-ik—l r- ' - 


- ArxkC^^-^BcdC=-^ 


k + l 


Here, we note that 
obtain 

(4.3) 




= At I = Ai t. By changing t and C each other, we 


pi 

ii,C -k BtC^-^Atx dr = 

Jo /c + i 


Differentiating the both sides with respect to leads to 
(4.4) 


B(;-k [ BtC^-^Atx dr = C^Bc^, 

Jo 


which is equivalent to (4.2). Note that (4.3) is equivalent to (4.4) due to (7o = 0. □ 
In the following subsections, we characterize the order conditions of energy¬ 
preserving CSRK methods in terms of the matrix M. 
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4.1. Energy-preserving integrators with = 1. First, we consider the 
simplest case = 1. 

Theorem 4.3. A CSRK method is energy-preserving and of order at least p = 2rj 
if the symmetric matrix M £ satisfies 

(4-5) ■■■ fe+s-i] M = ij, k = l,...,r], 

where all components of ik £ R* are zero except for the k-th component which is 1. 

Proof. We obtain = 1 by substituting k = 1 into (4.5), and this indicates that 
Cr = T. We consider an integrator satisfying the C{r]) condition. Since 



l/(s + fc - 1) 


this coincides with Cf/k = /k for any r if and only if 

-A- 1 1 M - i^ 

Ik fe+1 • ■ • fc+s-lj “ ‘fe • 


□ 

Remark 4. For the condition (4.5), it is possible to choose rj = s. In this case, 
(4.5) is eguivalent to 


1 1 
2 3 

1 _j_ 

_ S S+ 

and thus M can he written as M = H~^, where H is the Hilbert matrix defined by 
Hij = l/(z+j —1). In this case, the CSRK method coincides with the AVF collocation 
method of order 2s. 

4.2. Energy-preserving integrators with B(^ ^ 1. It is also possible to derive 
energy-preserving integrators which do not satisfy = 1. We here illustrate a 
derivation of 4-degree fourth order integrators. 

Let s = 4 and write [l \ 4 = v^. Then Bq and C^- are expressed as 


1 

s-l-l 


2s-l 


M = K 


Bc_ = 


c 

c' 



4 

r 

4 


V. 


If the C{2) assumption is satisfied, we have 

(«) = 
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This relation indicates that the polynomial degree of Ct is at most 2, because the 
polynomial degree in r of the left hand side is 4. We write 

= [r 2(1 — r) 0 O] , r ^ I 

among several possibilities for which d( = 1. It follows from (4.6) that 
^ ^ ^]^=[0 3r(l-r) 2(l-r)2]. 

Hence, if a symmetric matrix Af satisfies 


' 1 

1 

2 

1 

3 

1 

4 

M = 

V 2 ( 1 -r) 0 0 ’ 

r+2 

. 6 

r+3 

12 

r+4 

20 

r+5 
30 . 

0 3r(l — r) 2(1 — r)^ 


the method has order 4. As an example, if we select r = 0, this characterization 
reduces to 


'l 1 i i 


0200' 

2 3 4 

M = 


1111 
_3 4 5 6. 


0 0 0 2 


In this case, M can be expressed with three free parameters a, /3 and 7 : 


a 

1 

2 

/3 

1 ■ 

6 


'l 

1 

1 

1 

1 


1 

2 

4 

6 

8 

M = 

p 

1 

6 

7 

1 

10 


1 

1 

1 

1 

1 


1 

.6 

8 

10 

12. 



It seems that (a,/ 3 , 7 ) = (I,|;,i) are the simplest choices, and in this case M is 
calculated to be 


6 

5 

72 

5 

-36 

24 ' 

72 

5 

144 

5 

-48 

72 

-36 

-48 

720 

-720 

24 

72 

-720 

720 


As illustrated above, energy-preserving integrators can be obtained even for the 
cases Bq ^ 1. However, such cases require larger degrees, and thus are less practical 
than the cases Bc^ = 1. 

5. Parallel energy-preserving methods. In this section, we construct new 
energy-preserving methods which can be implemented more efficiently than the AVF 
collocation method. In the context of RK methods, it is known that the computational 
cost of solving implicit RK methods can be reduced if a RK matrix A has only real, 
distinct eigenvalues (see, e.g. [4, 13]). The key of the construction of new integrators is 
to apply similar idea to CSRK methods and put it together with the energy-preserving 
condition and order conditions. 
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5.1. Implementation of Runge Kutta methods with real eigenvalues. 

In general, in order to proceed with a RK process, we have to solve a system of 
nonlinear equations of size sN per each time step. Here we explain that if the matrix A 
has only real, distinct eigenvalues, the computational costs can be reduced, especially 
in a parallel architecture. 

We consider solving an s-stage implicit RK method 

S 

Y^=yo + h'^a^jf{Yj), i = 

1=1 

S 

yi = yo + h'^bJ{Yi). 

In order to compute the internal stages Yi,... ,Ys G we have to solve a system 
of nonlinear equations 

(5.1) ^{Y)=Y-es<S)yo-h{A^lN)F = 0 

where (8) denotes the Kronecker product. In G is an identity matrix and 



'Y{ 


■f 


'fXX 

Y = 


G R^"^, es = 


GR^ F = 



Xs. 


1 


fXs)_ 


The size of the system (5.1) is sN. We usually apply the simplified Newton method 
to solve the nonlinear system (5.1), and then obtain the iteration formula 

(5.2) {I,N-hA(^Jo)p‘ = -^Y^), Y‘+^=Y^ + p\ 1 = 0,1,2,..., 

where Jq denotes the Jacobian matrix, i.e. Jo = S'V^H{yo) for Hamiltonian systems. 

Below, we show that the linear system (5.2) of size sN can be computed efficiently 
if the matrix A has only real, distinct eigenvalues. 

If all eigenvalues of A are real and distinct, there exists a matrix T G such 

that 

T ^AT = diag(Ai,..., As), Ai,..., As G R. 

Let M = (IsN — hA 0 Jq) and define M hy M = {T~^ 0 In)M(T 0 In). Then it 
follows that 

M = diag(/Ar — hAi Jo ,... ,In — hAs Jo), 

where the right hand side expresses a block diagonal matrix. Therefore, the formula 

(5.2) is translated into 

YI+^=y^+p\ 1 = 0 , 1 , 2 ,..., 

where p^ is calculated based on the following relation 


(5.3) 

$(r') = (T-i®/^)$(y'), 

(5.4) 

II 

1 

(5.5) 

p' = {T-^0In)pK 
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The key point of this idea is that the linear system (5.4) of size sN consists of s linear 
systems of size N. 

Below we consider to what extent the computational cost is reduced. For both 
approaches, we equally pay operations to compute Jo- 

• Standard iteration based on (5.2). 

We compute the LU factors of {IsN — hA 0 Jq), where a small constant 
times (sN)^ operations are needed. Then we solve (5.2) iteratively until 
convergence, where 0{s^N'^) operations are needed per each iteration. 
Remark 5. Note that the east of the LU faetorisation can be further redueed 
by making use of complex conjugate eigenvalues pairs [13]. For example, the 
cost of the L U factorisation for the 2-stage Gauss method can be reduced from 
0{8N^) to OiAN^). 

• Iteration based on (5.4). 

We compute the LU factors of (/at — hAi Jq), ... ,{In — hX^Jo), where only 
0{sN^) operations are needed. They can be reduced to 0{N^) in a parallel 
architecture. Next we solve (5.4) iteratively until convergence, where only 
0{sN‘^) operations are needed, and they can be further reduced to 0{N‘^) 
if parallelism is available. Note that although we compute (5.3) and (5.5) 
additionally, the computational costs of these parts are just 0{N), and thus 
less important. 


5.2. Efficient continuous stage Runge Kutta methods and their imple¬ 
mentation. We apply the concept of the previous subsection for RK methods with 
real eigenvalues to CSRK methods. 

First, let us express W in (2.5) as 

S 

Yr = yolo{T) + ^ Ycih(T) 

2 = 1 

where ^^(t) is defined by 

s 

hir) = TT -—— (co = 0), f = 0, l,...,s. 

r\ • / • 

3^1 


We regard (2.5) as a system of nonlinear equations in terms of Wi i : Ws: 

(5.6) Yc,=yo-\-h[ Ac,,c/(lc) dC, i = l,...,s. 

Jo 

Here we merely evaluate (2.5) at the nodes t = ci,..., Cs. As this was done for RK 
methods, the expression (5.6) is further re-expressed in a form 

^'A,„c/(dc)dC 


(5.7) 


where 


$(y) = Y - es^yo - h 


/'A,,,c/(dl)dC 

■Jo 


= 0 , 


Yr, 


Yr. 


Y = 


e 
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Since the derivation of $ with respect to Y becomes 



• [\,,.,cJ{Y^)UC)dC 

Jo 

Jo 

/•I 

/■i 

/ A,,,cd(Fc)^i(C)dC 

LJo 

• / A,^,cJiyc)uc)dc 

Jo -1 


$'(y) = i,N-h 


the application of the simplified Newton method for solving (5.7) gives the iteration 
formula 

{I,N -hE0j)p^ = -^Y^), Y^+^=Y^+p\ 1 = 0 , 1 ,..., 

where the matrix i? e is dehned by 

(5.8) [ A,,,c^j(C)dC. 

Jo 

Therefore, if the matrix E oi a CSRK method has only real, distinct eigenvalues, the 
CSRK method can be implemented efficiently. The following theorem indicates that 
one does not have to be concerned about the Ci values when evaluating the eigenvalues 
ofE. 

Theorem 5.1. The eigenvalues of the matrix E defined in (5.8) are independent 
of the Ci values. 

Proof. We define a Vandermonde matrix 


V = 


Cl cf 
C2 (?2 


Note that the matrix V is non-singular, if the Ci values are distinct. 

Below, we show that V~^EV is independent of the Ci values. We express 
^r,c = Yli=i where Pi,..., are polynomials in C. Since Ac,x = c^PkiC), 

the matrix E is expressed as 

= /' E (c) dc = E (c) dc. 

0 7 . _1 7 _ 1 0 


k^l 




Therefore, it follows that 
(5.9) 


{V-^EVfi, = E dC ■ 4 = dC. 

fc=ido Jo 


In the second equality, the relation X]fe=i ^k^k{C,) = is used. □ 
The expression of V~^EV (5.9) can be further simplified. Since 


p^{0 = - [M^l 


M,, 
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where the matrix M is defined in (2.6), the component of 


V ^EV can be written as 


{V-^EV\, = P,(C)C^ dC = - [" [Mn 

Jo * Jo 

1 / ^ ^ M,s \ 

i\j + l j + sj' 

Hence, we have the following formula 


M,s] 


(^j + s-l 


(5.10) V-^EV = 

\ 

1 

2 

M 

1 1 

2 3 

1 1 

3 4 

1 

s+1 

1 

s-1-2 


1 

s_ 


1 1 
_ 1 s-\-2 

1 

2s 


dC 


In the following subsections, we construct efficient fourth and sixth order energy¬ 
preserving integrators by making use of the idea presented in this subsection. 


5.3. Three-degree fourth order integrators. 

5.3.1. Derivation. We derive fourth order energy-preserving integrators. Since 
the only fourth order CSRK method with degree two is the AVF method, we set the 
degree of CSRK methods to s = 3. We also assume that = 1. According to 
Theorem 4.3, a CSRK method is of order four if the symmetric matrix M satisfies 


1 ^ 

2 3 


1 

111 

M = 

1 

2 3 4 

i 3 “ 


1 


with a parameter a e M. When a ^ 7/36, M is analytically calculated to be 


ai -1- 4 

—6ai — 6 

6ai 

—6ai — 6 

36ai -f 12 

—36a 

6ai 

—36ai 

36ai 


M = 


The characteristic polynomial of the matrix (5.10): 


V-^EV = 


is found to be 


^(A) = - -A^ + 


1 

12 


(y.\ 


1 


ai = 


36a-7 


1 


111 

T 


2 3 4 

1 

M 

111 

2 


3 4 5 

1 


111 

3 


4 5 6 


A- 


ai 

600’ 


It follows that <f{X) = 0 has three real, distinct roots if 


> ^2^/3 -h —23/3 + -fv 0.7770503941. 

300 6 24 4 

Calculating the range of a is easy, and exactly the same evaluation was done by 
Butcher-Imran [4] in the context of symplectic methods. 
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5.3.2. Long time behaviour. Since the new method surely preserves the en¬ 
ergy, one may expect a good long-time behaviour. As a criterion of the long-time 
behaviour, we here consider the conjugate-symplecticity [11]. 

A one-step method of order p is said to be conjugate-symplectic up to order 
p -|- r (r > 0), if there exists such a change of coordinates z = xiu) that is 0{hP) close 
to the identity and '^h = X° o x~^ satisfies 

The modified equation of such a method is Hamiltonian up to the terms of size 
(^(/jP+r), |.]^g behaviour of the method is like that of a symplectic integrator 

on intervals of length 0{h~'^). See [11, Chapter VI.8] for more details. 

In [10, 14], it is proved that the AVF collocation method of order 4 (2s in gen¬ 
eral) is conjugate-symplectic up to order 6 (2s -I- 2), but not conjugate-symplectic 
up to a higher-order. It follows from Theorem 5.11 in Hairer-Zbinden [14] that the 
new method (of order 4) is conjugate-symplectic up to order 6 but not conjugate- 
symplectic up to a higher-order, independently of the choice of a. Therefore, one 
could expect that the new fourth-order method behaves like the fourth-order AVF 
collocation method in terms of conjugate-symplecticity. 

5.3.3. Efficiency. Although the computational cost in each time step of the de¬ 
rived integrators is reduced to about one-eighth (or one-fourth if the AVF collocation 
method is implemented by making use of complex conjugate eigenvalues pairs), the 
actual local truncation error of the new integrators is worse than the AVF collocation 
method. Therefore, we here carefully consider the efficiency of the proposed inte¬ 
grators by taking both truncation errors and computational costs into consideration. 
Below, we denote —ai/300 by 9. 

Table 1 shows coefficients of the elementary differentials with trees of order 5 
for the B-series expansions of the exact solution, the two-degree fourth-order AVF 
collocation method and the proposed method. It is estimated that when we use the 
same stepsize the error of the proposed method is approximately (600 -I- 1) times 
bigger than that of the fourth-order AVF collocation method. On the other hand, for 
large problems, the computation of the proposed method is 8 (or 4) times faster than 
the AVF collocation method per each time step. Thus it is fair to compare the local 
errors of the AVF collocation method with the stepsize 8h and the proposed method 
with stepsize h a.t t = to + 8h so that the overall computational costs are almost the 
same. The former is denoted by errAVF and the latter by errNew Then the ratio of 
the errors at t = to + 8h is roughly estimated to be 

8 • errNew _ 8 • (600 -k l)h^ _ 600 -k 1 
errAVF {8h)^ 4096 

Therefore, the proposed method is more efficient than the fourth-order AVF colloca¬ 
tion method (in the sense that the above ratio is less than 1) if 

273 

0.7770503941 < 0 < -= 68.25. 

4 

Note that if we consider the new method is only 4 times faster than the AVF collo¬ 
cation method, this range of 0 changes to 

0.7770503941 <9 <^ = 4.25, 

4 

and thus the proposed method is still efficient. 
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Table 1 

Coefficients of the elementary differentials with trees of order 5 for the B-series expansions of 
the exact solution, the fourth-order AVF collocation method and the proposed method. 


t 





y 

T 

f 

T 

1 

exact solution 

1 

5 

1 

10 

1 

15 

1 

30 

1 

20 

1 

20 

1 

40 

1 

60 

1 

120 

AVF collocation 

1 

T 

1 

TT7 

5 

TI 

5 

1 

1 

1 

1 

Tl 

1 

AVF collocation — exact 

0 

0 

1 

1 

0 

0 

0 

1 

1 

Tin 

Proposed method 

1 

Z 

1 

TTJ 

120-1-5 

— Tr~ 

120-1-5 

144 

1 

1 

1 

1-120 
— T2 — 

1-120 

144 

Proposed method — exact 

0 

0 

600-1-1 

yoo 

600-1-1 

Tib 

0 

0 

0 

600-1-1 

360 

600-t-l 

7‘lb 


5.4. Five-degree Sixth order integrators. We next consider the derivation 
of five-degree sixth order integrators. Since it turned out that intended integrators 
cannot be constructed with s = 4, we here set s = 5. We also assume that = 1. A 
CSRK method is of order at least six if the symmetric matrix M satisfies 


1 

1 

1 

1 

1 


1 


2 

3 

4 

5 


1 

1 

1 

1 

1 


1 

2 

3 

4 

5 

6 


1 

1 

1 

1 

1 

M = 

1 

3 

4 

5 

6 

7 

1 

4 

1 

5 

1 

6 

a 

P 


1 

1 

.5 

1 

6 

1 

7 

p 

7_ 


1 


with parameters a,/3,f € M. Recall that in the derivation of the three-degree fourth 
order integrators, it was easy to evaluate eigenvalues of the matrix (5.10), because 
there is only one free parameter. However, we find it increasingly difficult to do a 
similar estimate for sixth order integrators: the application of the Sturm theorem 
is too complicated due to three free parameters. Here, we only show that there are 
choices of parameters such that the matrix (5.10) has only real, distinct eigenvalues, 
by giving a concrete choice. 

It is difficult to express the matrix M analytically in terms of the parameters a, 
j3 and 7 . Thus, instead of requiring all eigenvalues of (5.10) to be real, we require real 
eigenvalues of 


We write 





a = 




Let 7 >(A) = det(A/ — N'). It is possible to find a parametric solution to the problem 
(/ 7 ( 0 ) = If)'{Q) = 0 in terms of a*, (5* and 7 *. One way of doing this is explained in 
Step 1 below. Here are three steps that need to be carried out to find a real root 
solution to ip. 
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1. We choose the following parameters 

, _ 1 + la? 3 + ac ^ 64 + 

“ “ 2800 ’ ^ ^ 4200 ’ ^ 44100 

where a and c are real or purely imaginary numbers, so that the first and 
second degree terms in (fi{X) are both 0 . 

2. It is found that 

y,(A) = A2(A3 + (d - 2)A2 - 6(5A + 125) 

where S = 18(21a — c)^ — 10. This polynomial has three real roots for 5 € 
[0,0.2144863035] where the upper bound is the real zero of 

5^ + 245^ + 1445 - 32 = 0. 

We choose the parameters a and c so that 5 £ [0,0.2144863035]. A possible 
choice is (a, c, 5) = (1, i). 

3. We add small perturbation to a or /I or 7 in the correct direction so that the 
two zero roots remain real. For example, reducing 7 by 10“^*^ gives the real 
roots 

-0.8831274189, -0.06591290801, 0.06591290335, 0.9185056830, 1.839542361 

Remark 6 . Although the five-degree sixth order integrator has been derived above, 
it is still diffieult to give a necessary and sufficient condition for the parameters a, fi 
and 7 such that the matrix (5.10) has only real eigenvalues. Therefore, we here leave 
the discussion about the efficiency untouched. 

6. Concluding remarks. In this paper, we have shown the energy-preserving 
condition of CSRK methods for Hamiltonian systems, characterized the order condi¬ 
tions and parallelizable condition in terms of the matrix M. As applications, we have 
derived new fourth and sixth order integrators which can be implemented more effi¬ 
ciently than the standard method, i.e. the AVF collocation method. In the derivation, 
the eigenvalues of the matrix E played an important role. 

We note several directions for future work. From a theoretical viewpoint, it is 
hoped that the conjecture discussed in Remark 3 will be solved. Considering other 
applications of the energy-preserving condition of CSRK methods would be interest¬ 
ing. We are currently trying to give a more systematic approach to deriving higher 
order efficient energy-preserving integrators. 
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